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ABSTRACT 

Zj \ It has recently been shown that the rate of angular momentum relaxation in 

nearly-Keplerian star clusters is greatly increased by a process termed resonant 
relaxation ( Rauch fc Trcmainc 1996]) ; it was also argued, via a series of scaling 
arguments, that tidal disruption of stars in galactic nuclei containing massive 
black holes could be noticeably enhanced by this process. We describe here the 
results of numerical simulations of resonant tidal disruption which quantitatively 
test the predictions made by Rauch & Tremaine. The simulation method is 
based on an A-body routine incorporating cloning of stars near the loss cone and 
a semi-relativistic symplectic integration scheme. Normalized disruption rates for 
resonant and non-resonant nuclei are derived at orbital energies both above and 
below the critical energy, and the corresponding angular momentum distribution 
functions are found. The black hole mass above which resonant tidal disruption 
is quenched by relativistic precession is determined. We also briefly describe the 



discovery of chaos in the Wisdom- Holman symplectic integrator applied to highly 
eccentric orbits and propose a modified integration scheme that remains robust 
under these conditions. 

?_i ■ We find that resonant disruption rates exceed their non-resonant counterparts 

by an amount consistent with the predictions; in particular, we estimate the net 
tidal disruption rate for a fully resonant cluster to be about twice that of its non- 
resonant counterpart. No significant enhancement in rates is observed outside 
the critical radius. Relativistic quenching of the effect is found to occur for hole 
masses M > Mq = (8 ± 3) x 10 7 Mq. The numerical results combined with the 
observed properties of galactic nuclei indicate that for most galaxies the resonant 
enhancement to tidal disruption rates will be very small. 

Subject headings: black hole physics — galaxies: active — galaxies: nuclei — 
stellar dynamics 
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1. Tidal Disruption in Galactic Nuclei 
1.1. Introduction 

It has long been speculated that massive black holes (MBHs) reside in the centers of 
many or most galaxies (e.g., [Lynden-Bell 19691) — a notion strongly supported by recent 
high-resolution observations of galactic nuclei QEckart fc Genzel 1996| ; |van der Marel et al 



1997| ; |Kormendy et al. 1996] ; |Kormendy fc Richstone 1995| and references therein); these 



observations also indicate that the density of stars rises continuously towards the center of 
the galaxy. It is nearly certain, therefore, that tidal disruption of stars by central MBHs 
actively occurs in many galaxies. 

Possible observable consequences of tidal disruption in galactic nuclei have been pro- 
posed by numerous authors. Hills (1975) first suggested that the ultimate power source 
for active galactic nuclei (AGNs) might be the accretion onto an MBH of gas from tidally 
disrupted stars (cf. the related model by [Stoeger, Pacholczyk, fc Stepinski 1992| ). Although 
subsequent calculations ( It 1 rank 1978| ; [McMillan, Lightman, fc Gohn 1981 ; Duncan fc Shapiro 



|1983| ) showed that this mechanism was incapable of sustaining typical AGN luminosities — 
the stellar densities and dispersions required implied that physical collisions would dominate 
the mass loss — tidal disruptions may nonetheless contribute significantly to the luminosity 
of some AGNs, particularly the less energetic variants such as Seyferts and LINERs (low 
ionization nuclear emission line regions). The X-ray outbursts observed in some Seyfert- 
like nuclei flGrupe et al. 19951 ; Bade, Komossa, fc Dahlem 1996|) , for example, could be the 



result of tidal disruptions. Additional observations that can be understood in terms of stel- 
lar disruptions include the variations in the strength and profile of Balmer lines seen in the 
Seyfert/LINER NGC 1097 (|Storchi-Bergmann et al. 1995|) and other active nuclei ([Eracleous 



|et al. 1995]) , and the absence of compact nuclear UV sources in ~ 80% of a sample of LIN- 
ERs (|Eracleous, Livio, fc Binette 1995| ). The latter, for instance, was interpreted using a 
'duty cycle' model in which the emission line region is powered by episodic accretion events 
(produced by tidal disruptions) lasting a few decades each but active only ~ 20% of the time 
on average. A model by Roos (1992) hypothesizes that broad-line clouds in AGNs are the 
remnants of tidal disruptions. It has also been suggested that the shell-like structure Sgr A 
East in the Galactic Center region may be the result of a disruption event occurring ~ 10 4 yr 
ago ( [Khokhlov fc Melia 1996 ). Other possible signatures of tidal disruptions include tran- 



sient flares in otherwise quiescent nuclei ( [Rees 1988D , high- velocity stars produced by tidal 



dissociation of binaries (the other star being left bound to the MBH; |Hills 198"8| ; [Hills 1991] ) , 
and nuclear enrichment in particular isotopic species due to the explosive nucleosynthesis 
which might occur in relativistic disruption events ( [Luminet fc Bar buy 1990] ). It has even 



been suggested (|Carter 1992|) that cosmological 7-ray bursts could be the result of extreme 
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disruption events in distant nuclei. 

Research into tidal disruption can be divided into two classes: studies of the dynamics of 
the tidally disrupted debris itself, and studies of the disruption-induced dynamical evolution 
of the nuclear star cluster. The former has been investigated by several groups (e.g., |Carter 



fe Luminet 1983| ; [Evans fc Kochanek 1989| ; |Laguna et al. 19931 ; [Kochanek 1994|) , their most 



basic finding being that approximately half of the disrupted debris remains bound to the hole 
(and is eventually accreted) with the remainder being ejected from the nucleus. In this paper 
we undertake an analysis of the latter type: an investigation of loss cone dynamics under the 
influence of resonant relaxation ( |Rauch fc Tremaine 1996| , hereafter RT96) , an enhancement 
to two-body relaxation that is especially pronounced in nearly-Keplerian stellar systems. 

Unless explicitly overridden, we will employ geometric units throughout this paper. In 
these units, G = c = M = 1 (where M is the hole mass) and the natural units of time, 
distance, and angular momentum are GM/c 3 , GM/c 2 , and GM/c, respectively. 



1.2. Review of Non- Resonant Loss Cone Dynamics 

The contribution of tidal disruptions driven by two-body relaxation to the dynamical 
evolution of star clusters has been examined in detail ( [Frank fc Rees 1976 ; Lightman fc 



Shapiro 1977| ; [Shapiro fc Marchant 1978| ; |Cohn fc Kulsrud 1978|) . The loss cone itself can 



be defined as the region of phase space for which the periapsis of the stellar orbit, r p , lies 
inside the star's tidal radius, r t « 2(M/m+) 1 / 3 R ic ; these orbits are normally highly radial. 
Restated in terms of angular momentum, loss cone orbits are those with L < L min , where 
L = [a(l — e 2 )] 1//2 is the angular momentum of a Keplerian orbit with semi-major axis a and 
eccentricity e, and L m i n = [(2 — r t /a)rt\ 1 / 2 rs (2r t ) 1 ^ 2 is the angular momentum of an orbit 
with periapsis r p = r t . It is useful to define a parameter q = (AL orb /L min ) 2 , where AL orb is 
the rms change in angular momentum (due to two-body relaxation) per orbital period; there 
are then two limiting cases to consider. For q <C 1, stars take many orbital periods to diffuse 
across the width of the loss cone and hence stars entering the loss cone disrupt within an 
orbital period; this is the "empty loss cone" or "diffusion" limit in which the phase space 
density drops exponentially (with scale length ~ AL orb ) in the region L < L min . For q ^> 1, 
stars can relax across the loss cone in less than an orbital period and thus only a fraction 
~ 1/q of the stars crossing the loss cone are disrupted (recall that L-space is effectively two 
dimensional in this problem) ; this is the "full loss cone" or "pinhole" limit for which the loss 
cone is dynamically unimportant and phase space densities inside and outside the loss cone 
are nearly equal. Since L min = L min (£) and AL orb m AL orb (£) (where £ = I /(2a) > is 
the orbital binding energy) — for realistic density profiles, the two-body relaxation time at 
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fixed £ depends only weakly on L — q g(£) and one can define a "critical energy," £ crit , for 
which q(£ CT it) = 1- Somewhat more loosely, one can define the boundary by a critical radius 

?"crit ~ 1/^crit- 

The qualitative influence of the loss cone on the overall structure of the cluster is now 
easy to see. For £ < £ crit (r ^ r crit ), q > 1 (the enclosed cluster mass fraction is larger, 
orbits are less nearly Keplerian, and hence they relax relatively faster) and the loss cone has 
minimal effect on the distribution function, f{£ , L). Outside the critical radius, therefore, the 
density profile will follow the classical p* oc r -7 / 4 profile ( [Bahcall fc Wolf 1976| ; this assumes 
an isotropic cluster and collisionless cusp growth — cf. IQuinlan, Hernquist, fc Sigurdssori 



1995[ ) and the cluster will remain locally isotropic, i.e., f(£,L) ~ /(£)• Conversely, for 
£ > £ crit (r ^ r crit ), q < 1 and the removal of stars by tidal disruption is significant; the 
density cusp begins to flatten inside r cr it, eventually turning over and then reaching zero 
at r t (the detailed calculations indicate that the flattening inside r cr i t begins logarithmically 
slowly — by r ~ 10r t , the cusp slope flattens by only ~ 0.2). In addition, phase space densities 
at fixed energy are severely depleted at low angular momenta: f(£, L) ~ for L < L min . For 
typical "cuspy" nuclear density profiles (p* oc r -7 with 1 < 7 < 2, which encloses the range 
observed in real galaxies; [Lauer et al. 1995Q , the total rate of disruptions is dominated by 
the contribution from stars with energies near £&&■ The simulations of Shapiro & Marchant 
(1978), for example, found the mean energy of disrupted stars to be ~ 3£ C rit (the median 
value, on the other hand, was near £ cr it)- 

An order of magnitude estimate for £ cr - lt follows easily from the definition g(£ C rit) = 
(AL or b/X m m) 2 = 1. First, let L max = a 1 / 2 denote the maximum angular momentum for an 
orbit with semi- major axis a. Then AL orb /L max ~ (torbAr) 1 ^ 2 and L min /L max ~ (2r t /a) 1 / 2 , 
where t T is the two-body relaxation time (e.g., |Spitzer fc Hart 1971| ). Assuming in addition 



R± oc mV 3 , p+ oc r 7 / 4 , and an orbit-averaged value for t r , we find 

8/27 

M 8 - 1/27 pf>, (1) 

where M 8 = M/(10 8 Mq) and p6 = p*(l pc)/(10 6 Mq pc -3 ). In physical units, the corre- 
sponding value of r crit = l/(2£ crit ) is 




-8/27 



r crit ~ 10 ^ j Mf / 27 pe 4/9 PC (2) 



We will parameterize the relative importance of tidal disruption in a cluster using a 
dimensionless function A (£), defined to be the fraction of stars at energy £ consumed in the 
loss cone each orbital period. Let Aq denote the theoretical value of A for a fully non-resonant 
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cluster. In the pinhole limit (g > 1), the loss cone is dynamically unimportant and Ao is 
simply the fractional area of the energy hypersurface occupied by the loss cone: 

Ao^f^V, q>l. (3) 



In the empty loss cone limit, q <C 1, the loss rate can be found by solving the Fokker-Planck 
equation (e.g., |(John fc Kulsrud 1978[ ), which yields 

a r 2 

UL Wh /L .y ?<1- (4) 

One result of our analysis will be a comparison of A(£) between equivalent resonant and 
non- resonant systems (§ |3.1|) . 



1.3. Resonant Relaxation Summary 

Under the usual two-body relaxation process, both the energy £ and angular momentum 
L of stellar orbits diffuse (i.e., random walk) slowly through phase space due to the mutual 
gravitational interactions with other stars in the cluster; on average, therefore, the variations 
in 8 and L, A£ and AL, grow with time oc t 1//2 , and the energy and angular momentum 
relaxation times coincide, tg ~ t^. As discussed in detail by RT96, resonant relaxation is 
a process by which angular momentum relaxation can be dramatically enhanced in clusters 
whose mean potential contains 'resonant' structure (energy relaxation is not affected). It 
operates most effectively in nearly-Keplerian star clusters, whose resonant structure consists 
of spherical symmetry and equality of the radial and azimuthal orbital frequencies. In this 
case, the changes in both the vector and scalar angular momenta grow approximately linearly 
with time — systematically faster than a random walk — and hence the angular momentum 
relaxation time is much shorter than the energy relaxation time: £l -C (possibly by several 
orders of magnitude). 

Because the nuclear cluster of every galaxy with a central MBH is nearly-Keplerian close 
to the black hole, resonant relaxation can have important consequences for the dynamics of 
galactic nuclei. Whenever ti, < tjj < tg, for instance (where t# is the Hubble time), resonant 
relaxation will act to isotropize the cluster in less than the age of the universe, even though 
the two-body relaxation time is much longer than tjj (as is typical for galactic nuclei). In 
the form of resonant dynamical friction, resonant relaxation can exert a controlling influence 
on the eccentricity evolution of a massive black hole binary (cf. IQuinlan fc Hernquist 1997| ); 
it can also rapidly erode the inclination of relatively massive bodies in nearly-Keplerian disk 
systems. For details, see RT96. 
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In this paper we examine in depth another consequence of resonant relaxation: resonant 
tidal disruption. In RT96 it was shown that the enhanced relaxation of L should lead to 
a large increase of the tidal disruption rate in the empty loss cone region (5 < 1), by a 
factor of roughly q~ 1 ^ 2 ] since disruption rates are not enhanced in the full loss cone region, 
however, they estimated that the overall disruption rate would increase only moderately (a 
factor of two or three at most). It was also argued that relativistic precession due to the 
MBH should effectively disable resonant tidal disruption for hole masses M^4x 10 7 Mq. 
(Resonant relaxation only operates on timescales shorter than the orbital precession time, 
which becomes dominated by relativistic effects as the tidal radius approaches the horizon.) 
However, resonant relaxation is a highly non-linear phenomenon, and it is difficult for analytic 
calculations to achieve better than order-of-magnitude accuracy; these numerical estimates 
are, therefore, rather crude. 

In the present work we discuss the results of detailed numerical simulations designed to 
quantitatively test these predictions. We calculate both the resonant enhancement to the 
tidal disruption rate and the associated change in the cluster distribution function. The 
influence of relativistic precession on the results is also investigated, yielding an improved 
estimate of the "quenching mass" above which resonant tidal disruption is disabled. We 
begin in § |2] with a description of the numerical method. Results are presented in § BL and 
concluding discussion is given in § [|. Appendices |A] and [FJ provide additional discussion of 
technical issues pertaining to the calculations. 



2. Simulation Method 

2.1. Numerical Strategy 

Numerical simulation of resonant tidal disruption is challenging for several reasons. 
First is the need to use N-body techniques in a case where the physical system under study 
contains many millions of stars. The non-resonant investigations described in § |1.2| avoided 
this problem by using a continuum approximation (the Fokker-Planck equation) to calculate 
the dynamics; in the resonant case, however, the Fokker-Planck equation is invalid — it is 
based on random-walk diffusion of L, whereas resonant relaxation leads to linear growth. 
Although it should be possible to modify the Fokker-Planck formalism to accurately account 
for resonant effects, this has not yet been accomplished; a particle-based approach is there- 
fore required. We ease the computational demands by differentiating stars into two types, 
test and background (§ |2.2|) , the former serving as the tracers of dynamic evolution and 
the latter representing the mass responsible for driving that evolution. Further reduction in 
computational effort was achieved using a specialized mass spectrum (and associated soft- 
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ening length) for the background stars, which allowed fewer stars to be included at large 
distances without significantly changing the physical relaxation time (§ |2.3|). 



A second constraint is the high dynamic range needed in L-space: determining the 
resonant modification to the cluster distribution function, f(L), requires adequately resolving 
/ on scales ~ L min ; in galactic nuclei, however, L max /L min ~ 10 5 for r ~ r crit . We deal with 
this problem in two ways: first, by restricting the simulations to a manageable subregion 
of L-space defined by L < L , where L min < L < L max ; and second, by cloning stars 
approaching L = L min to help reduce statistical fluctuations in this sparsely populated 
region of phase space (§ |2.4| ). Note that attempting to utilize a Fokker-Planck approach 
would also be complicated by these low phase space densities present near the loss cone; 
previous studies, by contrast, treated the loss cone as a boundary layer instead of resolving 
it explicitly. 

The radial dynamic range of the test orbits also places stringent demands on the choice 
of integrator — another consequence of the extreme mismatch between L min and L max , which 
implies that loss cone orbits are highly radial (in our simulations, e ^ 0.9999). As dis- 
cussed in § P2§, finding a suitable integration scheme proved an arduous task; although itself 



not trouble-free, the Wisdom-Holman symplectic mapping ( |Wisdom fc Holman 1991 ; this 



scheme is also known as the mixed-variable symplectic method) emerged as the clear win- 
ner. A slightly modified form of this method, which we term "semi-relativistic" symplectic 
integration, was used in simulations examining the influence of general relativity on resonant 
tidal disruption. 

A cornerstone of our analysis is the use of comparative simulations — i.e., simulations 
that were strictly identical apart from the strength of resonant effects — to determine the 
importance of resonant tidal disruption in a particular system. This 'differential measure- 
ment' approach removes as much as possible any systematic biases that may be introduced 
by the abstractions used to keep the problem computationally tractable. In any event, the 
relatively good agreement between the theoretical and measured values of quantities such as 
£crit (§ |3.1|) indicate that such systematic effects are probably small. 



2.2. Model Overview 

The basic computational model contained only two components: a central black hole 
and the surrounding nuclear cluster. Most simulations treated the black hole as a Newtonian 
point mass; those including relativistic corrections (§ |3.2|) assumed a Schwarzschild geometry. 
In both cases the black hole was held fixed at the center of the cluster — i.e., we neglect the 
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Brownian motion of the hole, which is dynamically unimportant in these systems (see § |4j). 
Black hole masses in different simulations rang ed from 10 6 to 1O 8 M (in the range of those 
inferred from observations) and remained fixed for the duration of the run — i.e., we neglect 
any evolution of the black hole due to accretion of tidal debris. 

The nuclear cluster was modeled as a halo of bound, 1 Mq stars with a distribution 
function of the form f(£) oc £ l l 4 for £ min < £ < £ max , and f(£) = elsewhere — in other 
words, as an isotropic cusp with a radial density profile oc r~ 7 / 4 in the range a min <^ r <^ 
a max , where a min = l/(2£ max ) and a max = l/(2£ min ) are the corresponding semi-major axes; 
this is the Bahcall & Wolf (1976) solution, valid for r ^> r t (see § |1.2| ). On physical grounds, 
Omin ^ r t an d a max ^ r h are required, where r h (~ 1/cr 2 for a cluster with velocity dispersion 
a) is the black hole's dynamical sphere of influence, outside of which the cusp solution need 
not apply. In practice, a min and a max were restricted to a narrower band centered on r crit 
in order to minimize the number of simulated stars; typical values were a m i n ~ r cr i t /30 
and a max ~ 10r cr i t . The absolute density was chosen to keep the mass fraction small at 
the critical radius, a requirement for resonant relaxation to be effective there; in most cases, 
p*(l pc) ~ 10 3 Mq pc -3 was used. Typical galactic nuclei have p*(l pc) ~ 10 4 — 10 6 Mq pc~ 3 
(cf- §§• 

Although the use of N-body techniques in our simulations was mandatory (§ |2.1|) , a 
fully self-consistent N-body integration would have been wasteful; since there is minimal 
dynamical evolution of the nucleus as a whole during the course of the simulations, following 
the detailed trajectory of every cluster star is unnecessary. For this reason the simulations 
separated the cluster into two types of stars, test and background, according to the treatment 
of their dynamical evolution. 

Conceptually the background stars represent the mass of the cluster, providing the 
perturbing force felt by the test stars; the test stars, by contrast, are the dynamic tracers 
of the cluster's equilibrium distribution and tidal disruption rate. Hence only the test stars' 
dynamical evolution is considered — they feel the forces of the black hole and the background 
stars (but not of the other test stars), and their orbits relax nearly as they would in a 
fully self-consistent integration; background stars feel only the central force and execute 
unperturbed Keplerian motion for the duration of the simulation. If there are n test stars 
and N ^> n background stars, the force calculation in this scheme is then 0(nN) instead of 
the much larger 0((n + iV) 2 ) for a self-consistent calculation. The detailed handling of each 
is explained in the following sections. 
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2.3. Background Cluster Specifics 



Given a min , a max , and p 6 = p*(l pc)/(10 6 Mq pc~ 3 ), the number of background stars, 
N, is known. Even at the rather low simulated densities, however, there were too many 
background stars to include individually — most clusters having N ~ 10 5 — 10 6 and the prac- 
tical limit being a few thousand. The simplest solution to this problem would be to increase 
the mass of every background star in proportion to their total overabundance; we adopted a 
more flexible strategy in which mass was made a function of the star's (constant) semi-major 
axis, the most massive "stars" occupying the largest, most distant orbits. This approach 
allowed the background cluster to retain both its fine-grained structure at small radii, where 
average interstellar distances are shortest, and its overall mass, which is dominated by the 
stars at large radii. The original mass density profile was preserved by concomitantly de- 
creasing the stellar number density at large radii. More precisely, the background stars were 
given a mass spectrum of the form 



In addition, for a > a c their number density was forced to be oc r -3 , resulting in a mass 
density oc r~ 7 / A at all radii, as desired; for a < a c , stars had their nominal mass and a 
number density oc r~ 7 / 4 . The initial orbital velocities of the more massive stars were also 
given a moderate tangential bias to prevent these masses from reaching small radii, where 
they would severely skew the local density; otherwise the orbits were drawn isotropically. 
The value of a c was determined by the requirement that iV masses, drawn at random from 
the above distribution, have (on average) a combined mass equal to that of the actual cluster 
(= m*N), where N was the number of background masses to include in the simulation; most 
runs used N = 2000. All simulations used = 1 Mq. 

Another consideration in setting up the background cluster was the amount by which to 
soften the point mass potentials of the background stars, which is necessary with constant 
timestep integrators (such as the Wisdom-Holman mapping) to maintain accuracy during 
close encounters (see § |2.4j below). Note first that using equal softening lengths for all stars is 
impractical — the background masses spanned several decades in radius, and using the same 
softening length for each would produce either negligible softening for the outermost masses 
or excessive softening for the innermost. We therefore decided to use a softening length 
proportional to the semi-major axis of the orbit. The more massive stars, however (those 
with a > a c ), were treated a bit differently in this regard. Because they actually represent 
a group of discrete stars which have been condensed into a single mass for computational 
efficiency (a type of "poor man's tree code"), the softening length for these objects should 
take account of their implied finite size; this in turn depends on their mass and average 




(5) 
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density. Based on the specific form of rh(a), the softening length b used for a particular 
background mass was: 



Thus b m 0.1a for a 3> a c , which provided a good fit to the actual mean potential of the 
massive stellar aggregates. In terms of b the softened potential was $ = — m/(r 2 + o 2 ) 1 / 2 . 

Finally, we comment on how the orbital integrations for the background stars were 
performed. As mentioned earlier the background stars follow strictly Keplerian orbits, and 
hence their orbital motion can be calculated analytically. If the integration timestep is small 
enough, however, direct integration will be faster since one analytic Kepler step is an order 
of magnitude slower than one direct step. To maximize efficiency, the code therefore used 
direct integration whenever a few such steps would yield sufficient accuracy; otherwise orbits 
were advanced analytically. Direct integrations employed a second-order leapfrog step which 
maintained a relative energy accuracy of 10 -3 at all times. 



ative performance of a wide variety of methods, including both direct integration schemes 
(such as leap frog, adaptive Runge-Kutta, and Bulirsch-Stoer) applied to either the ordinary 
or regularized equations of motion and specialized schemes (such as Ecke's method and the 
mixed-variable symplectic method) which model the motion as a perturbed Keplerian ellipse. 
The latter have the advantage of being able to use relatively large step sizes for particles 
undergoing nearly-Keplerian motion, as is true of the test stars. Symplectic integrators have 
the advantage that they are normally free of long-term growth of energy error, but only if a 
constant stepsize is used — hence they also have the disadvantage of being non-adaptive. As 
the extreme orbital eccentricities in this problem demand that the integrator be (in some 
sense) adaptive, non-adaptive direct methods such as leap frog or its higher-order variants 
were completely unusable. In the end none of the direct schemes proved to be competitive 
with the specialized ones. Of the two methods of the latter type, Encke's method turned 
out to be surprisingly inferior to the mixed-variable symplectic (hereafter, MVS) method, 
particularly since only the former can utilize an adaptive stepsize. Encke's method (e.g., 
Panby 1992|) works by recasting the equations of motion in terms of deviations from a fixed 
Keplerian reference orbit that is periodically reset (once per orbit, say). In theory, since 
these deviations are normally small and slowly- varying for nearly-Keplerian motion they can 
be integrated more easily than the original equations; in the limit of unperturbed Keplerian 




(6) 



2.4. Test Star Integrations 



Because of the special demands on the test star integrator 
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motion the deviations remain precisely zero and the method is exact for arbitrarily large 
stepsizes. The MVS method ( |Wisdom fc Holman 199T| ; cf. [Kinoshita, Yoshida, fc Nakai 



|1991[ ) can be thought of either as a symplectic analog of Encke's method in which the ref- 
erence orbit is reset every time step, or (better) as a specialization of the leap frog step in 
which the rectilinear motion between the impulsive velocity changes is replaced by a Keple- 
rian orbital segment; this method is also exact in the limit of unperturbed Keplerian motion. 
Our testing highlighted the fact that even if the absolute deviations integrated by Encke's 
method remain small, their derivatives can vary sharply — mainly near periapse. This im- 
plies that the deviation equations must be very carefully integrated when passing periapse, 
requiring the use of extremely small stepsizes for highly eccentric orbits; this probably ex- 
plains the method's poor performance in our problem. We therefore chose to use the MVS 
method for our calculations. Not even this scheme was problem-free, however, as it was 
found to inject a (remarkably!) large amount of spurious energy relaxation into the motion. 
Appendix [TJ discusses the origin of this behavior (which appears to be dynamical chaos in 
the MVS integrator) as well as the workaround we used to mitigate the problem. The final 
simulations were not adversely affected by the malady. A fringe benefit of the MVS scheme 
is the ease with which tidal disruptions can be detected; since the spatial motion between 
velocity perturbations is purely Keplerian, the minimum radius crossed during each step can 
be efficiently determined. 

As noted previously it was impractical to simulate all of L-space in our calculations — 
even resonant relaxation is too slow for this to be computationally feasible with N-body 
methods. We therefore focussed attention on a more manageable region of phase space 
around the loss cone, only considering test orbits with angular momenta L < Lq, where 
L m m <C L <C L max . Our estimates of the resonant modification to f(L) is thus confined 
to L < L ; since resonant effects are most dramatic near L min , however, this was not a 
significant restriction (see § |3|). There is, of course, nothing physically special about this 
boundary, and test stars inside the simulation boundary will freely relax across it unless 
manually restricted. We implemented this restriction by considering L = Lq to be a reflecting 
boundary — stars crossing outside L = L were simply removed from the simulation and 
replaced by clones (explained below) crossing the boundary towards the inside. In most 
simulations, L ^ 10 L min was attainable. 

To adequately resolve f(L) across the entire L < L subspace, the region was divided 
into a number of non-overlapping zones spaced logarithmically in L (typically ~ 20), each 
covering a factor of ~ v^2 in L. Determining f(L) (for a fixed energy E) was done by 
integrating a population of test stars, each with initial orbital energy E, until they achieved 
dynamic equilibrium; f(L) was then determined by computing the average occupancy of 
each zone (weighted by their relative areas). 
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The small relative areas enclosed by the innermost zones makes them subject to large 
statistical uncertainties (i.e., Poisson noise) due to the very small average occupancies within 
these zones. To achieve comparable relative errors in all zones, we introduced cloning of test 
stars into the simulations, which attempts to equalize the absolute occupancies between zones 
by duplicating stars entering the innermost zones. Conceptually our cloning procedure is 
very similar to that used by Shapiro & Marchant (1978) in their Monte Carlo loss cone 
simulations; it differs in that whereas their cloning occurred in energy space (allowing /(£) 
to be determined well above the critical energy, £ cr it), ours is done in L-space. The cloning 
process proceeded as follows. At the start of a simulation a number of zones were designated 
'n-tuple cloning zones,' meaning that any test star relaxing into such a zone from above (i.e., 
from an adjacent orbit with larger L) was replicated n times; subsequently, the n clones 
were integrated precisely like non-clones, except that any clone passing outside the zone in 
which it was created was removed from the simulation. Clones were, however, free to relax 
inwards (i.e., towards more radial orbits with smaller L) — possibly creating clones of their 
own, and so on; those suffering tidal disruption were also removed without being replaced. 
Disrupted non-clones, on the other hand, were replaced by new test stars lying just inside 
the reflecting boundary. Note that although the number of clones is dynamic, the number of 
non-clones (and thus the average test star "potential" ) remains fixed. The system is therefore 
dynamically closed and the zonal occupations come into statistical equilibrium after a few 
crossing times. (Recall that since test stars do not feel the gravitational forces of the other 
test stars, cloning introduces no dynamical feedback into the system.) We found that great 
care was needed in creating the clones to avoid biasing the calculated distribution function. 
In particular, by construction all stars undergoing cloning have dL/dt < (their orbits 
became more radial during the preceding integration step), and not only must the newly- 
created clones share this property, the distribution of their dL/dt values must match that of 
the uncloned stars. In practice this was achieved by creating candidate clones, integrating 
them backwards one timestep, and rejecting any that did not cross the cloning boundary 
during this step (as the original star had); otherwise the only physical variables the clones 
had in common with their parents were their energies £ and (scalar) angular momenta L. 

We tested the reflection and cloning schemes by performing a preliminary simulation in 
which tidal disruption was artificially turned off (by simply not checking for it). Since the 
perturbing potential provided by the background cluster is nearly isotropic, in equilibrium 
the test stars should exhibit a flat distribution function, f(L) = const. The calculated f(L) 
was constant to within « 3%, which is less than the statistical uncertainty of the results 
presented in § ||. Other checks — e.g., that the mean density profile of the background cluster 
was indeed oc r~ 7//4 — were also done and all returned nominal results. 
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3. Results 

The consequences of resonant tidal disruption were examined through a series of com- 
parative resonant and non-resonant simulations that were otherwise identical. The non- 
resonant results corresponding to a particular resonant simulation were created by redoing 
the resonant run while artificially adding orbital precession (amounting to ~ 1 rad per orbit) 
to the test stars' motion to force the precession time down to only a few orbital periods, 
which is short enough to extinguish resonant effects; operationally this was done by adding 
a small amount of precession each integration step. This technique allowed a "differential 
measurement" of the disruption rates and equilibrium distributions to be made, which is 
more robust and reliable than attempting to calibrate and match absolute results obtained 
from only qualitatively-similar runs. 

3.1. Newtonian Simulations 

This set of simulations was designed to explore the gross differences in disruption rates 
and dynamical equilibrium between resonant and non-resonant clusters. The runs used a hole 
mass M = 10 6 Mq , a value appropriate for small galactic nuclei such as the Galactic Center 
or M32 (cf. Table ^j. In this case, r t ~ 100 is well outside the horizon and relativistic effects 
are unimportant; hence these simulations were performed within an entirely Newtonian 
framework using the standard MVS algorithm. 

The main simulations consisted of four sets of comparative runs, each at a fixed orbital 
energy; three of the four were at an energy above the critical energy, where resonant en- 
hancements should be prominent. Phase space parameters for these calculations are shown 
in Table [TJ. The table lists the orbital energy of the test stars, £ , the maximum angular 
momentum followed in the simulation, Lq (i.e., the location of the reflecting boundary — see 
§ |2.4j) , the maximum possible angular momentum at the given energy, L max , the average 
amount of momentum relaxation suffered per orbit, AL or b, and the typical relaxation suf- 
fered by resonantly relaxed stars over a full precession time, AL prec ; subscripts "res" and 
"nr" refer to the values for resonant and non-resonant runs, respectively. Also given are 
the results for the normalized tidal disruption rates, A, with estimated 1 — a errors given 
in parentheses (Ao is the 'baseline' expectation; see § |1.2| ). The background cluster in these 
simulations had parameters (§ |2.3| ) N = 2000, p$ = 10 -3 (p* oc r~ 7 / 4 ), £ min = 2 x 10~ 9 , and 
£max = 5 x 10~ 7 ; it was also found empirically to have £ crit ~ 2 x 10~ 8 , in good agreement 
with the estimate given in eq. ([!]). 

The strong resonant enhancement to the local disruption rate is clearly evident in Ta- 
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ble [T]. To within the accuracy of the results (~ 25% at 1 — cr), the ratio of the resonant to 
non-resonant tidal disruption rate is given simply by 

Ares J -^min/^-^orb; AZ/ or b — L m [ ni , . 

\ ~ II A T > T ^ ' 

This agrees qualitatively with the estimate given by RT96, and quantitatively restricts the 
previously undetermined numerical coefficient to be very close to unity. Combining this result 
with the analytic estimate for A nr (§ |1.2| ) implies that the overall rate of tidal disruptions in 
a resonant nucleus will be at most three times that for an equivalent non-resonant nucleus 
(note: for p+ oc r~ 7 / 4 , it is easily shown that AL or b is roughly proportional to S' 1 ), and 
probably only a factor of two greater. The precise value depends on the density profile, 
the strength of resonant effects near £ CI i t , the inner radius at which disruptions significantly 
deplete the cusp, and the radius at which relativistic effects become relevant. 

A comparison of the equilibrium angular momentum distribution functions is shown in 
Figure ^ each panel corresponding to a specific value of £ (or AL or b, marked); all curves 
have been normalized to have unit area in the range [0,L max ]. For each curve, the errorbar 
for the outermost point includes the overall uncertainty in normalization (resulting from the 
need to extrapolate the curves from Lq to L max ) as well as the statistical uncertainty in the 
average bin occupancy, whereas the remaining points include only the latter. Uncertainty in 
the normalization of the resonant curves was estimated by splitting the difference between 
two extreme assumptions for how to extrapolate f(L), namely: (1) f(L) = const. = f(L ), 
or (2) df /d(\nL) = const. = df /d(\nL) nT for L > L . Note that in agreement with previous 
studies, the non-resonant f(L) are well fit by a straight line in semi- log coordinates (outside 
the loss cone); for these curves the extrapolation error was taken to be the formal uncertainty 
in the corresponding linear least-squares fit. The strong resonant increase in phase space 
density near the loss cone (which disappears outside the critical radius) is again clearly 
visible. 

Figure § compares the resonant (solid) and non-resonant (dotted) f{L) for a run with 
L ~ L max ; a linear least-squares fit to the latter (dashed line) is also plotted. The figure 
suggests that for large L the two curves are essentially equal — the resonant enhancement 
becoming significant only near the loss cone. This behavior can be understood qualitatively 
as follows. Recalling that even resonant relaxation proceeds as a random walk on timescales 
longer than the orbital precession time (§ |1.3| ), one would expect relaxation over distances 
much larger than AL prec — the effective mean free path in L-space — to proceed at almost 
identical rates in both resonant and non-resonant systems; their distribution functions should 
therefore be nearly coincident in this region. Near the loss cone, on the other hand (i.e., for 
L <^ AL prec ), the effective diffusion in the resonant cluster is much higher; here, therefore, the 
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resonant f(L) will be flatter than its non- resonant counterpart, and hence will increasingly 
dominate over the latter as the loss cone is approached. Inside the loss cone both curves 
fall exponentially with a scale length of AL or b. This line of reasoning explains the gross 
features visible in Figure |2|; in particular, note the apparent turnover in / res for L ^ 3L min 
(this particular simulation had AL 

prec ^ 2_ZL mm ). 



3.2. Relativistic Simulations 

We examined the importance of relativistic precession through two series of simula- 
tions. In each series, the value of AL OT ^/ L min was held fixed while the black hole mass 
was varied between relatively small values, for which relativistic effects were negligible, and 
larger masses for which they dominated the precession. The relativistic simulations were 
performed in a conceptually simple manner — by merely replacing the Kepler stepper used 
in the usual MVS scheme (§ |2.4|) with one integrating motion in the Schwarzschild (non- 
rotating black hole) spacetime. This is equivalent to including the relativistic corrections 
to the Newtonian black hole (i.e., point mass) potential in the "unperturbed motion" part 
of the symplectic integrator instead of in the perturbation potential. Although straightfor- 
ward, this approach does not treat time dilation effects completely accurately — time is still 
considered an absolute, universal quantity. Since such effects are noticeable only during the 
extremely brief passage through periapse, however, they can be safely neglected in our case. 
(In fact, our method is probably overkill; artificially adding precession at the relativistic rate 
would likely have been sufficient.) Because of this comingling of formalisms, we refer to this 
scheme as "semi-relativistic" symplectic integration; for details on implementation of the 
Schwarzschild stepper (which incorporates a high-order approximation to geodesic motion 
in the Schwarzschild metric) see Appendix 

Figure |3] displays the disruption rate A as a function of black hole mass for test stars 
undergoing strongly resonant relaxation (AL orb /L min = 0.025) in the non-relativistic case. 
The points A(M — > 0) and A(M — > oo) correspond to purely Newtonian resonant and non- 
resonant simulations, respectively. The figure shows very clearly the quenching effect of 
relativistic precession at large hole masses: whereas the scale-free Newtonian simulations 
predict a factor of ~ 30 enhancement in the tidal disruption rate due to resonant effects (cf. 
Table [I]), for a physical hole mass of (say) 10 6 Mq, the relativistic simulations show only 
a factor of ~ 5 net increase. The "quenching mass" at this energy is thus approximately 
10 6 Mq. Near the critical energy, on the other hand, relativistic effects from a 10 6 Mq 
black hole would not noticeably wash out the (much smaller, but more robust) resonant 
enhancement to A; accomplishing this requires a much shorter precession period — and hence 
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Fig. 1. — A comparison at several energies of the normalized L-space distribution func- 
tions, including approximate 1-a error bars, for equivalent resonant (solid curves with square 
points) and non-resonant (dotted curves with triangular points) clusters. The large, outer- 
most error bars include the uncertainty in overall normalization; the remainder show only the 
formal statistical uncertainties in each point. In the empty loss cone limit, AL or b/-L m m *C 1, 
resonant phase space densities are much higher near the loss cone, implying a large increase 
in the local tidal disruption rate relative to non-resonant clusters (cf. Table [1]). 
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Fig. 2.— A comparison of resonant (solid curve with square points) and non-resonant 
(dotted curve with triangular points) distribution functions near L max (marked). The dashed 
line through the non- resonant f(L) is a linear (in semi-log coordinates) least-squares fit to 
the data. For L 3> AL prec ~ 2, the two curves appear to merge; only nearer the loss cone is 
the resonant enhancement to f(L) significant. 
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a larger black hole. 

The second set of relativistic simulations were designed to estimate the hole mass above 
which all resonant tidal effects disappear. Three distinct groups of runs were performed, 
all with a fixed value of AL or b/L m in = 0.3: (1) non-relativistic, non-resonant; (2) non- 
relativistic, resonant; and (3) relativistic, resonant. We will define the quenching mass 
as the black hole mass for which half the total resonant enhancement to the overall tidal 
disruption is relativistically suppressed. In this case, the estimate is made more difficult by 
the relatively small resonant excess near the critical energy as well as the large array of runs 
required. From the results of ~ 20 runs spanning hole masses of 10 7 to 10 8 , our best estimate 
for the quenching mass Mq is 

M Q = (8 ±3) x 1O 7 M . (8) 
This is also close to the rough estimate of 4 x 10 7 Mq given in RT96. 



4. Discussion 

We have shown that resonant tidal disruption strongly enhances the disruption rate in 
the empty loss cone limit, by a factor of ~ L min /AL; the total disruption rate for the entire 
cluster is approximately double the non- resonant rate. We also found that for hole masses 
M ^ 8 x 1O 7 M , the latter increase is eliminated by relativistic precession; sufficiently far 
into the diffusion limit, this quenching effect is also present for smaller black holes. These 
results quantitatively confirm the order-of-magnitude estimates given in RT96 and highlight 
the potential importance of resonant relaxation in determining the structure of galactic nuclei 
at small radii. 

Other processes capable of modifying loss cone dynamics include black hole wander 
and tidal capture. Black hole wander ( |Lin fc Tremaine 1980| ; | Young 1977] ; |Bahcall fc Wolt 



|1976| ) — the result of two-body relaxation acting on the hole itself — is a Brownian motion-like 
process which causes the hole to undergo minor excursions from the true dynamical center 
of the nucleus. Like resonant tidal disruption, BH wander is never important in the full 
loss cone regime — not even the existence of a loss cone, much less its precise location, is of 
consequence in this case. In the empty loss cone case, BH wander can modestly enhance 
the local non- resonant disruption rate by a logarithmic factor ( [Young 1977| ) ; however, we do 
not expect resonant relaxation and BH wander to significantly reinforce each other in this 
regard. This is because resonant nuclei have relatively low densities and hence large r cr i t ; 
the low velocity dispersions here then imply loss cone impact parameters b t ~ L min /t> ^> r w , 
where r w is the rms BH wandering radius (even in the relatively dense Galactic Center, for 



-19- 



AL/L . =0.025 

' mm 



C 

X 




10 



log(M/M ) 



Fig. 3. — A plot of the fractional disruption rate per orbital period, A, as a function of the 
black hole mass, M, for a series of otherwise identical resonant clusters. For very massive 
holes, the resonant enhancement to A is quenched — the result of an increasing amount of 
relativistic precession near the loss cone as the tidal radius moves closer to the horizon. When 
relativistic precession is negligible (leftmost point), the disruption rate is ~ 30 times higher 
than in the corresponding non- resonant cluster (rightmost point). In this particular case, 
half of this increase is gone by M ~ 10 6 Mq; near the critical energy (where AL orb /L min ~ 1) 
the quenching mass is larger, ~ 8 x 10 7 Mq (see § |3.2|) . 
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example, b t ~ 3r w near r crit ). At smaller radii, by contrast, the tightly bound cusp stars 
have orbital periods short compared to the wandering timescale, and therefore follow the 
hole adiabatically. There could be an intermediate region in which BH wander and resonant 
effects reinforce each other, but in general this does not happen. 

Tidal capture, on the other hand, will further enhance resonant disruption rates just 
as it enhances non-resonant tidal disruption. Tidal capture (piener et al. 1995 ; |Novikov 



Pethick, fc Polnarev 19921 |Rees 1988| |Frank fc Rees 1976]) occurs when a star passes close to 



but outside of, its tidal radius; although it escapes immediate destruction, the tides raised 
on the star dissipate enough orbital kinetic energy inside the star to substantially reduce the 
semi-major axis of the orbit (the pericentric distance remains nearly fixed). This process 
continues until either the orbit circularizes or the star disrupts; in the case of a massive 
black hole, the studies indicate that the latter will occur. (We remark that the consequence 
of tidal capture in £-space is conceptually similar to that of resonant relaxation in L-space: 
a systematic drift that overwhelms the diffusion induced by two-body effects.) Thus tidal 
capture effectively increases the value of L min (by a factor of w 1.5; Novikov et al. 1992); 
the result is a further doubling of the stellar disruption rate. Since this gas is, on average, 
tightly bound to the hole, the increase in accretable gas will be even greater. This further 
improves the ability of the loss cone to power at least low-level activity in galactic nuclei. 

How important a role will resonant tidal disruption play in normal galactic nuclei? 
Table |2| lists the nuclear properties of several black hole candidates, including estimates of the 
critical radius, r crit , and the cluster mass fraction at that location, /i cr it. The table suggests 
that typical galactic nuclei are non-resonant at their critical radii (recall that resonant effects 
dominate only for mass fractions <^ 0.1); resonant enhancement to the overall disruption 
rates in these nuclei will therefore be quite small. At radii a few times smaller than r crit , 
however, most galaxies in the table are dominated by resonant relaxation, and these regions 
can be resolved in several instances. Because erosion of the density cusp inside r crit progresses 
so slowly, however (see § |1.2|), its detection is unlikely even after accounting for resonant 



effects. We therefore conclude that resonant tidal disruption will not visibly alter nuclear 
density profiles at currently accessible angular resolutions. 

Another astrophysical system in which resonant loss cone effects are important is the 
Oort cloud. Although several processes contribute, the dynamical evolution of Oort cloud 
comets in the outer Solar System (r ~ 10 4 AU) is controlled mainly by the Galactic tidal 
field ( |Wiegert fc Tremaine 19971 [Heisler fc Tremaine 1986] ). The tidal perturbation to the 
nearly-Keplerian potential of the Solar System resonantly relaxes the angular momentum of 
the comet orbits, occasionally removing so much of it that the comets plunge into the inner 
Solar System, where they may become visible as long period comets. The loss cone in this 
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situation arises from the ability of the giant planets (particularly Jupiter and Saturn) to 
scatter incoming comets onto either hyperbolic or relatively tightly bound orbits — in either 
case causing the comets to be "lost" from the Oort cloud. As for nuclear clusters around 
MBHs, the rate at which comets reach the scattering loss cone — as well as the equilibrium 
Oort cloud distribution function at small L (which directly determines the frequency of 'new' 
long period comets) — is set by resonant loss cone dynamics. 

The relative difficulty of the numerical calculations highlights the generic advantage of 
using a continuum approach in this class of problems, and indicates that the development of 
a resonant Fokker-Planck formalism certainly merits further attention. The inherently non- 
linear nature of resonant relaxation and the associated complexity of its analytic description, 
however, suggests that it is likely to be a highly non-trivial undertaking! Further exploration 
of this issue is beyond the scope of this paper. 

We thank Scott Tremaine for helpful discussions and a careful reading of the manuscript. 
This work was supported by NSERC and a Jeffrey L. Bishop Fellowship to K. R. 



A. Semi-Relativistic Symplectic Integration 



As discussed in § |3.2| , some simulations incorporated relativistic corrections to the test 
stars' orbital motion to determine the hole masses for which resonant tidal disruption is 
shut off by relativistic precession. Operationally this was accomplished by replacing the 
Kepler stepper normally used in the MVS integration scheme, which solves for motion in 
the Newtonian one-body problem, with a Schwarzschild stepper approximating the motion 
around a Schwarzschild (i.e., non-rotating) black hole. In this appendix we describe the 
construction of the Schwarzschild stepper itself. 

In the Newtonian case all five orbital elements — the semi-major axis, a, eccentricity, 
e, inclination, i, argument of periapse, u, and argument of the ascending node, T — are 
constants of the motion; hence orbits are closed and do not precess in the absence of per- 
turbations. In the Schwarzschild case, continued spherical symmetry (which ensures that 
the orbital motion is planar) combined with energy and (scalar) angular momentum con- 
servation imply that four of these elements — a, e, i, and T — remain constant and well- 
defined even for highly relativistic orbits; to, however, is not fixed but precesses by an 
amount Au k, 67r/[a(l — e 2 )] rad per orbit in the weak field limit. (The exact expres- 
sion is Aco = 2[(3K((3 2 5e) - tt], where 5 = l/[o(l - e 2 )], (3 = 2/[l - 26(3 - e)] 1 / 2 , and K(m) 
is the complete elliptic integral of the first kind.) 
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Table 1. Resonant and Non- Resonant Tidal Disruption Rates 



^xrit 


T a 

-^0,nr 


r a 

-^0,res 


T a 

^max 


AL or b a 


Af a 

^-^prcc 


A x 10 7 


A nr x 10 7 


Arcs X 10 7 


8 


2.8 


8 


125 


0.022 


~ 3 


0.06 


0.18(0.04) 


5.5(1.3) 


3 


8 


23 


205 


0.14 


~ 6 


0.8 


1.4(0.2) 


9(1.5) 


1 


23 


32 


354 


0.80 


~ 12 


~ 25 


18(2) 


35(6) 


0.3 


64 


91 


614 


5.5 


~ 24 


26 


26(1.1) 


26(1.8) 



a In units of L mi 



Table 2. Loss Cone Properties of Nearby Black Hole Candidates 



Galaxy 


D (Mpc) 


M 8 


P6 


n a 


Tent (") 


^crit (PC) 


pcrit 


Milky Way 


0.01 


0.025 


1 


1.8 


7 


0.3 


0.3 


M32 


0.7 


0.02 


0.2 


2 


0.1 


0.4 


0.3 


M31 


0.7 


0.3 


0.1 


1 


1.5 


5 


0.4 


NGC 3377 


9.9 


0.8 


0.2 


1.3 


0.25 


12 


2 


NGC 3115 


8.4 


20 


0.6 


1.75 


8 


350 


4 


M87 


15.3 


30 


0.01 


1.2 


8 


650 


4 



p*(r) = 10 6 pe(r/l pc)-"M pc~ 3 
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A Kepler stepper advances the orbital position x and velocity v by a time At through 
solution of some form of Kepler's Equation, which determines the orbital phase as a function 
of time. The standard form of Kepler's Equation for elliptical orbits is 



M = E - e sin E, 



(Al) 



where M = 2n(t — £ p )/£ rb is the mean anomaly, E is the eccentric anomaly (related to the 
radius by r = a(l — ecosE)), t p is a time of pericentric passage, and t OT ^ = 2ixc?l 2 is the 
orbital period. Conceptually, a relativistic geodesic stepper advances the 4-position X and 
4-momentum P over some interval of coordinate time At; however, since the Schwarzschild 
metric is static it is sufficient in this case to use just the 3 spatial components of position and 
momentum (or velocity) to advance orbits, even ones passing near the horizon. One practical 
disadvantage of using the true relativistic 3- velocity is that calculation of the orbital elements 
requires setting up and solving a cubic polynomial, which is time-consuming (as it must be 
redone every time step). For this reason we created a 'semi-relativistic' Schwarzschild stepper 
which returns the true orbital 3-position but an effective 3-velocity for which application of 
the simpler Newtonian formulae yield the correct (relativistic) orbital elements; in the weak 
field limit, both of these vectors become identical to their non- relativistic counterparts x 
and v. Our semi-relativistic stepper can thus be used as an exact replacement for the Kepler 
stepper in the integration of moderately-relativistic orbits. 

Implementing a Schwarzschild stepper requires an analog of Kepler's Equation which is 
valid for relativistic orbits; generalized expressions for t or b and the true anomaly, v (defined 
as the physical angle traversed in the orbital plane since periapse), are also needed. Rauch 
(1997) has derived these formulae for orbits in the Kerr (rotating black hole) spacetime 
accurate to second order in S (see above), assuming 5 <C 1. For Schwarzschild the results 
(with errors 0(5 3 )) simplify to 



f orl) = ■>-« '' 2 {\ ■ ( - 



1 + 35 + 



54 + 3e z 



X+ 1 + 



9 + —e cos y 
4 



(A2) 

5\Sesmx, (A3) 



and M = E - esin£ + AE, where AE = f (x - E)(l - e 2 f' 2 5 2 and 
e = e jl - (£) 1 + (-10 + 18e 2 + 15v / T^) j 



(A4) 



Here x is the "relativistic anomaly," defined according to r = a(l — e 2 )/(l + ecosx) — hence 
X advances by exactly 2n from one periapse to the next, regardless of the amount of orbital 
precession (which is Au = u(2n) — 2n ps 6n5, in agreement with the value given above); 
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for strictly Keplerian orbits (i.e., 8 — > 0), x — v - Note that since < e < e < 1, the 
modified Kepler's Equation can be solved economically by finding the solution for AE — > 
and then correcting for the finite value of AE. We found the resulting stepper to be only a 
factor of two slower than its Newtonian counterpart and to give excellent results (compared 
to accurate integrations of the Schwarzschild problem) for radii r ^ 10, at which point the 
systematic error in the precession rate was ~ 25%. 

B. Chaos in MVS Integration of Nearly- Radial Orbits 

By construction the MVS integration scheme (see § |2.4| ) is exact in the limit of unper- 
turbed Keplerian motion. As a symplectic mapping, the method might also be expected to 
be free of systematic, long-term growth of energy error in the slightly perturbed problem 
(i.e., for nearly-Keplerian motion) — this property being one of the key advantages of such 
algorithms. Applied to our simulations, however, the MVS integrator was found to exhibit 
unbounded error growth even for very small perturbations, in some cases causing test star 
orbits to become unbound in as few as 100 orbital periods (the physical energy relaxation 
timescale being orders of magnitude longer). As we explain below, the instability is related 
to the very high eccentricities of the test star orbits; qualitatively similar integrations (e.g., 
tracing the orbital evolution of Oort cloud comets) are therefore likely to suffer from the same 
problem. In this section we briefly discuss the problem, and propose a modified integration 
scheme which remains robust under such circumstances. 

Our integrations are somewhat unique for two reasons: the test star orbits are all highly 
eccentric (generally, e ^ 0.999), and the perturbations from the background cluster are 
very nearly constant over scales of order the pericentric distance — the perturbations are 
spatially well-resolved with a stepsize At/t m b ~ 10~ 3 , taking the form of a small, fixed tidal 
force near periapse. Hence there is no a priori reason to integrate using a stepsize small 
enough to resolve periapse (At/t orb ~ (1 — e) -3 / 2 ~ 10~ 5 ), which would be impractical 
in any case — recall that symplectic integrators do not easily admit variable stepsizes; this 
is (in principle!) the great advantage of the MVS method for nearly-Keplerian problems. 
For stepsizes not resolving periapse, however, it was found that significant jumps in the 
energy, centered around periapse, occurred each orbit. This produced a linearly growing 
energy error — even in simulations where the background masses were held fixed, for which 
energy should be strictly conserved. In our case, it was possible to work around the problem 
by artificially truncating the cluster at a relatively large radius — ~ a/30 for test stars 
with semi-major axis a — for which the residual tidal force at periapse was small enough 
that physical energy relaxation dominated that spuriously produced by the MVS method; 
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angular momentum relaxation is nearly unchanged as it is controlled by the largest-scale 
fluctuations in the background cluster. Although this proved adequate, it would be more 
elegant (and generically more useful) to modify the MVS scheme so that it did not suffer 
from this problem; although a detailed treatment is beyond the scope of this paper, we will 
mention one possible solution. 

Since the perturbing force is almost constant near periapse — where effectively all of 
the erroneous relaxation occurred — an excellent model analogy to our situation is provided 
by the Stark problem, in which the perturbing force is strictly constant in magnitude and 
direction over the entire orbit. Besides its simplicity, the Stark problem has the added 
advantages of being integrable and having a conserved energy; this guarantees that any 
energy relaxation or chaotic behavior found in numerical integrations of the problem are 
an artifact of the integrator itself. Sample integrations using the MVS scheme confirmed 
that for stepsizes resolving periapse, there was no systematic growth in the energy error, as 
nominally expected from a symplectic scheme; unresolved orbits, on the other hand, generally 
had erratic, unbounded growth of energy error, even for extremely small perturbations. 
More direct (though informal) evidence for dynamic chaos in the integrator was found in 
the sensitivity of the relaxation curves to the precise initial conditions used, which was not 
present in the stable integrations. 

The close correspondence between the Stark problem and our simulations suggests that 
a more robust integrator might result by replacing the Kepler stepper used in the MVS 
method with a Stark stepper solving that problem exactly (which can be done in parabolic 
coordinates using elliptic functions and integrals). In contrast to the standard MVS scheme, 
which is exact in the limit of zero perturbations, the Stark integrator is exact in the limit of 
a constant (possibly large) perturbation force. A brief investigation revealed that the new 
method, as expected, does not suffer from the aberrant relaxation present in MVS scheme- 
even for large stepsizes that do not resolve the passage through periapse. On the other hand, 
the practicality of a Stark integrator depends on the speed of the Stark stepper itself, which 
will be much slower (and more tedious to implement) than the Kepler stepper it replaces. 
Based on the encouraging initial tests, an efficient implementation of the proposed Stark 
integrator is currently being undertaken; detailed results will be reported elsewhere. 
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